home *** CD-ROM | disk | FTP | other *** search
/ The Arsenal Files 4 / The Arsenal Files 4 (Arsenal Computer).ISO / ham / sattrk31.tgz / sattrack-3.1.tar / SatTrack / src / sattrack / satcalc.c < prev    next >
C/C++ Source or Header  |  1995-03-16  |  29KB  |  846 lines
  1. /******************************************************************************/
  2. /*                                                                            */
  3. /*  Title       : satcalc.c                                                   */
  4. /*  Author      : Manfred Bester                                              */
  5. /*  Date        : 03Mar92                                                     */
  6. /*  Last change : 15Mar95                                                     */
  7. /*                                                                            */
  8. /*  Synopsis    : Auxiliary routines for the satellite tracking program       */
  9. /*                'sattrack'.                                                 */
  10. /*                                                                            */
  11. /*                                                                            */
  12. /*  SatTrack is Copyright (c) 1992, 1993, 1994, 1995 by Manfred Bester.       */
  13. /*  All Rights Reserved.                                                      */
  14. /*                                                                            */
  15. /*  Permission to use, copy, and distribute SatTrack and its documentation    */
  16. /*  in its entirety for educational, research and non-profit purposes,        */
  17. /*  without fee, and without a written agreement is hereby granted, provided  */
  18. /*  that the above copyright notice and the following three paragraphs appear */
  19. /*  in all copies. SatTrack may be modified for personal purposes, but        */
  20. /*  modified versions may NOT be distributed without prior consent of the     */
  21. /*  author.                                                                   */
  22. /*                                                                            */
  23. /*  Permission to incorporate this software into commercial products may be   */
  24. /*  obtained from the author, Dr. Manfred Bester, 1636 M. L. King Jr. Way,    */
  25. /*  Berkeley, CA 94709, USA. Note that distributing SatTrack 'bundled' in     */
  26. /*  with ANY product is considered to be a 'commercial purpose'.              */
  27. /*                                                                            */
  28. /*  IN NO EVENT SHALL THE AUTHOR BE LIABLE TO ANY PARTY FOR DIRECT, INDIRECT, */
  29. /*  SPECIAL, INCIDENTAL, OR CONSEQUENTIAL DAMAGES ARISING OUT OF THE USE OF   */
  30. /*  THIS SOFTWARE AND ITS DOCUMENTATION, EVEN IF THE AUTHOR HAS BEEN ADVISED  */
  31. /*  OF THE POSSIBILITY OF SUCH DAMAGE.                                        */
  32. /*                                                                            */
  33. /*  THE AUTHOR SPECIFICALLY DISCLAIMS ANY WARRANTIES, INCLUDING, BUT NOT      */
  34. /*  LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A   */
  35. /*  PARTICULAR PURPOSE. THE SOFTWARE PROVIDED HEREUNDER IS ON AN "AS IS"      */
  36. /*  BASIS, AND THE AUTHOR HAS NO OBLIGATIONS TO PROVIDE MAINTENANCE, SUPPORT, */
  37. /*  UPDATES, ENHANCEMENTS, OR MODIFICATIONS.                                  */
  38. /*                                                                            */
  39. /******************************************************************************/
  40.  
  41. #include <stdio.h>
  42. #include <math.h>
  43.  
  44. #ifndef STDLIB
  45. #include <stdlib.h>
  46. #endif
  47.  
  48. #include "satglobalsx.h"
  49. #include "sattrack.h"
  50.  
  51. /******************************************************************************/
  52. /*                                                                            */
  53. /* getTrueAnomaly: gets mean and true anomaly                                 */
  54. /*                                                                            */
  55. /******************************************************************************/
  56.  
  57. void getTrueAnomaly(tmArg)
  58.  
  59. double tmArg;
  60.  
  61. {
  62.     getOrbitNumber(tmArg);                          /* also gets mean anomaly */
  63.     getSatPrecession();
  64.     kepler();
  65.  
  66.     return;
  67. }
  68.  
  69. /******************************************************************************/
  70. /*                                                                            */
  71. /* getOrbitNumber: gets the orbit number and the mean anomaly                 */
  72. /*                 both orbit number and mean anomaly count from the perigee  */
  73. /*                                                                            */
  74. /******************************************************************************/
  75.  
  76. void getOrbitNumber(tArg)
  77.  
  78. double tArg;
  79.  
  80. {
  81.     double dT, dT2, orbFract;
  82.  
  83.     dT          = tArg - epochDay;
  84.     dT2         = dT * dT;
  85.     curMotion   = epochMeanMotion + 2.0*decayRate*dT + 6.0*decayRateDot*dT2;
  86.  
  87.     refOrbit    = (double) epochOrbitNum + epochMeanAnomaly / TWOPI; /* [rev] */
  88.     curOrbit    = refOrbit + curMotion * dT;
  89.  
  90.     orbitNum    = (long) curOrbit + 1L;
  91.     orbFract    = modf(curOrbit,&dummyd);
  92.     orbitFract  = orbFract * 100.0;
  93.     meanAnomaly = orbFract * TWOPI;
  94.  
  95.     return;
  96. }
  97.  
  98. /******************************************************************************/
  99. /*                                                                            */
  100. /* getSatPrecession: calculates precession of the satellite's orbit           */
  101. /*                                                                            */
  102. /*                   for reference see:                                       */
  103. /*                   "The Satellite Experimenter's Handbook", M. Davidoff,    */
  104. /*                   The ARRL, Newington                                      */
  105. /*                                                                            */
  106. /******************************************************************************/
  107.  
  108. void getSatPrecession()
  109.  
  110. {
  111.     double satPrecFact;
  112.  
  113.     getSemiMajorAxis();
  114.     checkPerigeeHeight();
  115.  
  116.     satPrecFact = pow((EARTHRADIUS / semiMajorAxis), 3.5) /
  117.                           SQR(1.0 - SQR(eccentricity)) * CDR;
  118.  
  119.     raanPrec    = -9.95 * satPrecFact * cosInclination;
  120.     perigeePrec =  4.97 * satPrecFact * (5.0 * SQR(cosInclination) - 1.0);
  121.  
  122.     return;
  123. }
  124.  
  125. /******************************************************************************/
  126. /*                                                                            */
  127. /* kepler: solves Kepler's equation iteratively                               */
  128. /*                                                                            */
  129. /******************************************************************************/
  130.  
  131. void kepler()
  132.  
  133. {
  134.     register double eccen, eccAnom, meanAnom, diffAnom;
  135.  
  136.     eccen    = eccentricity;
  137.     meanAnom = meanAnomaly;
  138.     eccAnom  = meanAnomaly;                              /* initial guess     */
  139.  
  140.     do
  141.     {
  142.         diffAnom = (eccAnom - eccen * sin(eccAnom) - 
  143.                     meanAnom) / (1.0 - eccen * cos(eccAnom));
  144.         eccAnom -= diffAnom;
  145.     }
  146.     while (fabs(diffAnom) >= CASR);                     /* 1 arcsec precision */
  147.  
  148.     if (fabs(eccAnom - PI) < CASR)
  149.         trueAnomaly = PI;
  150.     else
  151.         trueAnomaly = 2.0 * atan(sqrt((1.0 + eccen) / (1.0 - eccen)) * 
  152.                              tan(eccAnom / 2.0));
  153.  
  154.     trueAnomaly = reduce(trueAnomaly,ZERO,TWOPI);
  155.  
  156.     return;
  157. }
  158.  
  159. /******************************************************************************/
  160. /*                                                                            */
  161. /* getSubSatPoint: calculates sub-satellite point in geodetic coordinates     */
  162. /*                                                                            */
  163. /******************************************************************************/
  164.  
  165. void getSubSatPoint(sspObject)
  166.  
  167. int sspObject;
  168.  
  169. {
  170.     double geocenPosSSP[3], geocenPosEquat[3];
  171.     double rObj, rEquat, lng, rEff, a, c, h, e2, z, lat, lat1;
  172.     int    i;
  173.  
  174.     if (sspObject == SAT)
  175.     {
  176.         for (i = 0; i <= 2; i++)
  177.             geocenPosSSP[i] = satPosGl[i];
  178.     }
  179.  
  180.     if (sspObject == SUN)
  181.     {
  182.         for (i = 0; i <= 2; i++)
  183.             geocenPosSSP[i] = sunPosGl[i];
  184.     }
  185.  
  186.     rObj = absol(geocenPosSSP);
  187.  
  188.     geocenPosEquat[0] = geocenPosSSP[0];
  189.     geocenPosEquat[1] = geocenPosSSP[1];
  190.     geocenPosEquat[2] = ZERO;
  191.  
  192.     rEquat = absol(geocenPosEquat);
  193.  
  194.     geocenPosEquat[0] /= rEquat;
  195.     geocenPosEquat[1] /= rEquat;
  196.  
  197.     lng = gasTime - atan2(geocenPosEquat[1],geocenPosEquat[0]);
  198.     lng = reduce(lng,-PI,PI);
  199.  
  200.     z   = geocenPosSSP[2];                                            /* [km] */
  201.     lat = atan(z/rEquat);
  202.     e2  = 2.0 * EARTHFLAT - SQR(EARTHFLAT);
  203.     a   = EARTHRADIUS;                                                /* [km] */
  204.  
  205.     do
  206.     {
  207.         lat1 = lat;
  208.         c    = 1.0 / sqrt(1.0 - e2 * SQR(sin(lat1)));
  209.         lat  = atan((z + a * c * e2 * sin(lat1))/rEquat);
  210.       
  211.     } while (fabs(lat1 - lat) > ONEPPM);
  212.  
  213.     if (sspObject == SAT)
  214.     {
  215.         h            = rEquat / rObj * a / cos(lat) - a * c;
  216.         rEff         = a + h;
  217.  
  218.         satLat       = lat;
  219.         satLong      = lng;
  220.  
  221.         satHeight    = rObj- rEff;
  222.         satCrashFlag = (satHeight < 0.0) ? TRUE : FALSE;
  223.     }
  224.  
  225.     if (sspObject == SUN)
  226.     {
  227.         sunLat       = lat;
  228.         sunLong      = lng;
  229.     }
  230.  
  231.     return;
  232. }
  233.  
  234. /******************************************************************************/
  235. /*                                                                            */
  236. /* getStateVector: calculates satellite's state vector, i.e. position and     */
  237. /*                 velocity in the mean equatorial coordinate system (M50)    */
  238. /*                                                                            */
  239. /*                 This function also provides essential calculations for the */
  240. /*                 orbit number and squint angle calculations and therefore   */
  241. /*                 needs to be called even when the state vector has been     */
  242. /*                 calculated with the SGP4 model already.                    */
  243. /*                                                                            */
  244. /******************************************************************************/
  245.  
  246. void getStateVector(timeArg)
  247.  
  248. double timeArg;
  249.  
  250. {
  251.     double pVec[3], qVec[3];
  252.     double vxOrb, vyOrb, vzOrb;
  253.     double xObs, yObs, zObs, xObsOrb, yObsOrb, zObsOrb;
  254.     double tmp, cosArgPerigee, sinArgPerigee, cosRaan, sinRaan;
  255.     double phi, cosPhi, sinPhi;
  256.     double xOrb = ZERO;
  257.     double yOrb = ZERO;
  258.     double zOrb = ZERO;
  259.     int    i;
  260.  
  261.     if (propModelType == TLEMEAN)
  262.     {
  263.         getTrueAnomaly(timeArg);
  264.         curArgPerigee  = epochArgPerigee + (timeArg - epochDay) * perigeePrec;
  265.         curRaan        = epochRaan + (timeArg - epochDay) * raanPrec;
  266.         curMeanMotion  = curMotion;
  267.         curInclination = inclination;
  268.     }
  269.  
  270.     else
  271.     {
  272.         meanAnomaly    = meanAnomalyX;
  273.         trueAnomaly    = trueAnomalyX;
  274.         curArgPerigee  = curArgPerigeeX;
  275.         curRaan        = curRaanX;
  276.         curMeanMotion  = curMotionX;
  277.         curInclination = curInclinationX;
  278.     }
  279.  
  280.     cosTrueAnomaly = cos(trueAnomaly);
  281.     sinTrueAnomaly = sin(trueAnomaly);
  282.  
  283.     cosArgPerigee  = cos(curArgPerigee);
  284.     sinArgPerigee  = sin(curArgPerigee);
  285.  
  286.     cosRaan        = cos(curRaan);
  287.     sinRaan        = sin(curRaan);
  288.  
  289.     satRadius      = semiMajorAxis * (1.0 - SQR(eccentricity))
  290.                      / (1.0 + eccentricity * cosTrueAnomaly);
  291.  
  292.     /* the next lines are needed in any case for the squint angle calculation */
  293.  
  294.     if (propModelType == TLEMEAN || attitudeFlag)
  295.     {
  296.         xOrb    =  satRadius * cosTrueAnomaly;   /* position in orbital plane */
  297.         yOrb    =  satRadius * sinTrueAnomaly;
  298.         zOrb    =  ZERO;
  299.     }
  300.  
  301.     /* calculate state vector if not yet calculated from SGP4/SDP4 model */
  302.  
  303.     if (propModelType == TLEMEAN)
  304.     {
  305.         tmp     =  sqrt(GMSGP / (semiMajorAxis * (1.0 - SQR(eccentricity))));
  306.  
  307.         vxOrb   = -tmp * sinTrueAnomaly;         /* velocity in orbital plane */
  308.         vyOrb   =  tmp * (cosTrueAnomaly + eccentricity);
  309.         vzOrb   =  ZERO;
  310.  
  311.         pVec[0] =  cosArgPerigee*cosRaan - sinArgPerigee*sinRaan*cosInclination;
  312.         pVec[1] =  cosArgPerigee*sinRaan + sinArgPerigee*cosRaan*cosInclination;
  313.         pVec[2] =  sinArgPerigee*sinInclination;
  314.  
  315.         qVec[0] = -sinArgPerigee*cosRaan - cosArgPerigee*sinRaan*cosInclination;
  316.         qVec[1] = -sinArgPerigee*sinRaan + cosArgPerigee*cosRaan*cosInclination;
  317.         qVec[2] =  cosArgPerigee*sinInclination;
  318.  
  319.         for (i = 0; i <= 2; i++)
  320.         {
  321.             satPosGl[i] =  xOrb * pVec[i] +  yOrb * qVec[i];
  322.             satVelGl[i] = vxOrb * pVec[i] + vyOrb * qVec[i];
  323.         }
  324.     }
  325.  
  326.     /* calculate squint vector if satellite's attitude is specified */
  327.  
  328.     if (attitudeFlag)
  329.     {
  330.         phi     = reduce(lasTime - curRaan - HALFPI,ZERO,TWOPI);
  331.         cosPhi  = cos(phi);
  332.         sinPhi  = sin(phi);
  333.  
  334.         xObs    = siteVecGl[0];
  335.         yObs    = siteVecGl[1];
  336.         zObs    = siteVecGl[2];
  337.  
  338.         xObsOrb = (-sinArgPerigee*cosInclination*cosPhi + 
  339.                 cosArgPerigee*sinPhi) * xObs + 
  340.               ( sinArgPerigee*cosInclination*sinPhi + 
  341.                 cosArgPerigee*cosPhi) * yObs -
  342.                 sinArgPerigee*sinInclination * zObs;
  343.  
  344.         yObsOrb = (-cosArgPerigee*cosInclination*cosPhi - 
  345.                 sinArgPerigee*sinPhi) * xObs + 
  346.               ( cosArgPerigee*cosInclination*sinPhi - 
  347.                 sinArgPerigee*cosPhi) * yObs - 
  348.                 cosArgPerigee*sinInclination * zObs;
  349.  
  350.         zObsOrb =  -sinInclination*cosPhi * xObs +
  351.                     sinInclination*sinPhi * yObs + 
  352.                 cosInclination * zObs;
  353.  
  354.         satSquintGl[0] = -xOrb - xObsOrb;
  355.         satSquintGl[1] = -yOrb - yObsOrb;
  356.         satSquintGl[2] = -zOrb - zObsOrb;
  357.     }
  358.  
  359.     return;
  360. }
  361.  
  362. /******************************************************************************/
  363. /*                                                                            */
  364. /* getSiteVector: computes the site position and velocity in the mean         */
  365. /*                equatorial coordinate system. The matrix 'siteRotMat' is    */
  366. /*                used by getTopoVector() to convert geocentric coordinates   */
  367. /*                to topocentric (observer-centered) coordinates.             */
  368. /*                                                                            */
  369. /******************************************************************************/
  370.  
  371. void getSiteVector()
  372.  
  373. {
  374.     double rotVec0[3], rotVec1[3], rotVec2[3];
  375.     double lat, cosSiteLat, sinSiteLat, siteRA, cosSiteRA, sinSiteRA;
  376.     double sqf, a, c, s, h, vecLen;
  377.     int    i;
  378.  
  379.     lat         = siteLat;
  380.     cosSiteLat  = cos(lat);
  381.     sinSiteLat  = sin(lat);
  382.  
  383.     siteRA      = lasTime;
  384.     cosSiteRA   = cos(siteRA);
  385.     sinSiteRA   = sin(siteRA);
  386.  
  387.     sqf = SQR(1.0 - EARTHFLAT);
  388.  
  389.     a = EARTHRADIUS;
  390.     c = 1.0 / sqrt(SQR(cosSiteLat) + sqf * SQR(sinSiteLat));
  391.     s = sqf * c;
  392.     h = siteAlt;
  393.  
  394.     /* this vector is the geocentric position of the site [km] */
  395.  
  396.     sitePosGl[0] = (a*c + h) * cosSiteLat * cosSiteRA;
  397.     sitePosGl[1] = (a*c + h) * cosSiteLat * sinSiteRA;
  398.     sitePosGl[2] = (a*s + h) * sinSiteLat;
  399.  
  400.     /* this vector is not affected by the non-spherical Earth */
  401.  
  402.     siteVelGl[0] = -SIDRATE * sitePosGl[1];
  403.     siteVelGl[1] =  SIDRATE * sitePosGl[0];
  404.     siteVelGl[2] =  ZERO;
  405.  
  406.     /* set up the three unit vectors for the rotation matrix */
  407.  
  408.     rotVec1[0] = -sinSiteRA;
  409.     rotVec1[1] =  cosSiteRA;
  410.     rotVec1[2] =  ZERO;
  411.  
  412.     vecLen = absol(sitePosGl);
  413.  
  414.     for (i = 0; i <= 2; i++)
  415.         rotVec2[i] = sitePosGl[i] / vecLen;
  416.  
  417.     cross(rotVec2,rotVec1,rotVec0,-1);
  418.  
  419.     /* this matrix is the geocentric site rotation matrix */
  420.  
  421.     siteRotMatGl[0][0] = rotVec0[0];
  422.     siteRotMatGl[0][1] = rotVec0[1];
  423.     siteRotMatGl[0][2] = rotVec0[2];
  424.     siteRotMatGl[1][0] = rotVec1[0];
  425.     siteRotMatGl[1][1] = rotVec1[1];
  426.     siteRotMatGl[1][2] = rotVec1[2];
  427.     siteRotMatGl[2][0] = rotVec2[0];
  428.     siteRotMatGl[2][1] = rotVec2[1];
  429.     siteRotMatGl[2][2] = rotVec2[2];
  430.  
  431.     return;
  432. }
  433.  
  434. /******************************************************************************/
  435. /*                                                                            */
  436. /* getRange: calculates satellite range                                       */
  437. /*                                                                            */
  438. /******************************************************************************/
  439.  
  440. void getRange()
  441.  
  442. {
  443.     double satDist[3];
  444.     int    i;
  445.  
  446.     for (i = 0; i <= 2; i++)
  447.         satDist[i] = satPosGl[i] - sitePosGl[i];
  448.  
  449.     satRange    = absol(satDist);
  450.  
  451.     rangeRate   = ((satVelGl[0] - siteVelGl[0]) * satDist[0] + 
  452.                    (satVelGl[1] - siteVelGl[1]) * satDist[1] +
  453.                    (satVelGl[2] - siteVelGl[2]) * satDist[2]) / satRange;
  454.  
  455.     satVelocity = absol(satVelGl);
  456.  
  457.     return;
  458. }
  459.  
  460. /******************************************************************************/
  461. /*                                                                            */
  462. /* getTopoVector: converts from geocentric mean equatorial coordinates to     */
  463. /*                topocentric coordinates                                     */
  464. /*                                                                            */
  465. /*                topocentric coordinate system used: right-handed            */
  466. /*                +X = south,  +Y = east,  +Z = up                            */
  467. /*                                                                            */
  468. /******************************************************************************/
  469.  
  470. void getTopoVector(objGeo,pobjTopo)
  471.  
  472. double objGeo[3], (*pobjTopo)[3];
  473.  
  474. {
  475.     double objDistT[3], objTopo[3];
  476.     double objDistNorm;
  477.     int    i;
  478.  
  479.     for (i = 0; i <= 2; i++)
  480.         objDistT[i] = objGeo[i] - sitePosGl[i];
  481.  
  482.     objDistNorm = absol(objDistT);
  483.  
  484.     for (i = 0; i <= 2; i++)
  485.         objDistT[i] /= objDistNorm;
  486.  
  487.     multMatVec(objDistT,objTopo,siteRotMatGl);
  488.  
  489.     for (i = 0; i <= 2; i++)
  490.         (*pobjTopo)[i] = objTopo[i];
  491.  
  492.     return;
  493. }
  494.  
  495. /******************************************************************************/
  496. /*                                                                            */
  497. /* getSunPhaseAngle: calculates phase angle between the ground station and    */
  498. /*                   the Sun, as seen from the satellite                      */
  499. /*                                                                            */
  500. /******************************************************************************/
  501.  
  502. void getSunPhaseAngle()
  503.  
  504. {
  505.     double satDistGnd[3], satDistSun[3], satNormGnd, satNormSun, arg;
  506.     int    i;
  507.  
  508.     for (i = 0; i <= 2; i++)
  509.     {
  510.         satDistGnd[i] = sitePosGl[i] - satPosGl[i];
  511.         satDistSun[i] = sunPosGl[i]  - satPosGl[i];
  512.     }
  513.  
  514.     satNormGnd    = absol(satDistGnd);
  515.     satNormSun    = absol(satDistSun);
  516.     arg           = scalar(satDistGnd,satDistSun) / satNormGnd / satNormSun;
  517.     sunPhaseAngle = fabs(acos(arg));
  518.  
  519.     return;
  520. }
  521.  
  522. /******************************************************************************/
  523. /*                                                                            */
  524. /* getAziElev: calculates object's vector position in the local geodetic      */
  525. /*             frame, as well as azimuth and elevation angles                 */
  526. /*                                                                            */
  527. /******************************************************************************/
  528.  
  529. void getAziElev(switchFlag)
  530.  
  531. int switchFlag;
  532.  
  533. {
  534.     double  localVec[3], aziVec[3], aziPerp[3], refrVec[3];
  535.     double  vectLen, azimuth, elevation, arg;
  536.     double  refIndex, refCorr, diffCorr, airMass;
  537.     int     i;
  538.  
  539.     if (switchFlag == SAT)
  540.     {
  541. /*  if astronomical precession is needed the following three lines have to be */
  542. /*  activated and the fourth line needs to be commented out                   */
  543. /*
  544.         getPrecMatrix();
  545.         multMatVec(satPosGl,satPosPrec,precMatrix);
  546.         getTopoVector(satPosPrec,localVec);
  547. */
  548.  
  549.         getTopoVector(satPosGl,localVec);
  550.     }
  551.  
  552. /*  for the Sun and the Moon the apparent place is calculated already and     */
  553. /*  astronomical precession is not needed                                     */
  554.  
  555.     if (switchFlag == SUN)
  556.         getTopoVector(sunPosGl,localVec);
  557.  
  558.     if (switchFlag == MOON)
  559.         getTopoVector(moonPosGl,localVec);
  560.  
  561.     aziVec[0]  = localVec[0];
  562.     aziVec[1]  = localVec[1];
  563.     aziVec[2]  = ZERO;
  564.  
  565.     vectLen    = absol(aziVec);
  566.  
  567.     aziVec[0] /= vectLen;
  568.     aziVec[1] /= vectLen;
  569.  
  570.     if (aziVec[0] == ZERO)
  571.         azimuth = (aziVec[1] > ZERO) ? HALFPI : THREEHALFPI;
  572.     else
  573.         azimuth = PI - atan2(aziVec[1],aziVec[0]);
  574.  
  575.     if (azimuth < ZERO)
  576.         azimuth += PI;
  577.  
  578.     arg = localVec[2];
  579.  
  580.     if (arg >  1.0) arg =  1.0;      /* protect against out-of-range problems */
  581.     if (arg < -1.0) arg = -1.0;
  582.  
  583.     elevation = asin(arg);
  584.  
  585.     /* apply correction for atmospheric refraction for visible tracking */
  586.  
  587.     getRefraction(atmPressure,ambTemperature,relHumidity,
  588.                   WAVEL,REFWAVEL,EARTHRADIUS,SCALEHEIGHT,
  589.                   siteAlt,elevation,&refIndex,&refCorr,&diffCorr,&airMass);
  590.  
  591.     aziPerp[0] =  aziVec[1];
  592.     aziPerp[1] = -aziVec[0];
  593.     aziPerp[2] = ZERO;
  594.  
  595.     cross(aziPerp,localVec,refrVec,1);
  596.  
  597.     vectLen = refCorr;
  598.  
  599.     for (i = 0; i <=2; i++)
  600.         localVec[i] += vectLen * refrVec[i];
  601.  
  602.     vectLen = absol(localVec);
  603.  
  604.     for (i = 0; i <=2; i++)
  605.         localVec[i] /= vectLen;
  606.  
  607.     arg = localVec[2];
  608.  
  609.     if (arg >  1.0) arg =  1.0;      /* protect against out-of-range problems */
  610.     if (arg < -1.0) arg = -1.0;
  611.  
  612.     elevation = asin(arg);
  613.  
  614.     if (switchFlag == SAT)
  615.     {
  616.         for (i = 0; i <= 2; i++)
  617.             localVecSatGl[i] = localVec[i];
  618.  
  619.         satAzimuth   = azimuth;
  620.         satElevation = elevation;
  621.     }
  622.  
  623.     if (switchFlag == SUN)
  624.     {
  625.         for (i = 0; i <= 2; i++)
  626.             localVecSunGl[i] = localVec[i];
  627.  
  628.         sunAzimuth   = azimuth;
  629.         sunElevation = elevation;
  630.     }
  631.  
  632.     if (switchFlag == MOON)
  633.     {
  634.         for (i = 0; i <= 2; i++)
  635.             localVecMoonGl[i] = localVec[i];
  636.  
  637.         moonAzimuth   = azimuth;
  638.         moonElevation = elevation;
  639.     }
  640.  
  641.     return;
  642. }
  643.  
  644. /******************************************************************************/
  645. /*                                                                            */
  646. /* satEclipse: calculates the eclipses of the satellite                       */
  647. /*                                                                            */
  648. /* this function assumes that the Sun is point-like rather than extended      */
  649. /* and that the Earth is perfectly spherical                                  */
  650. /*                                                                            */
  651. /******************************************************************************/
  652.  
  653. int satEclipse(elev)
  654.  
  655. double elev;
  656.  
  657. {
  658.     double satPara, satPerp, satDist, alpha, beta;
  659.  
  660.     satDist = absol(satPosGl);
  661.  
  662.     satPara = scalar(satPosGl,sunPosGl) / sunDist;
  663.     satPerp = sqrt(SQR(satDist) - SQR(satPara));
  664.  
  665.     /* case 1: satellite is on side of the Earth facing the Sun */
  666.  
  667.     if (satPara > ZERO)   /* satellite is always in daylight */
  668.     {
  669.         if (elev < ZERO)
  670.             return(DAY);
  671.         
  672.         if (elev >= ZERO && sunElevation <= TWILIGHT)
  673.             return(VISIBLE);
  674.         else
  675.             return(DAY);
  676.     }
  677.  
  678.     alpha = asin(EARTHRADIUS / sunDist);
  679.     beta  = atan(satPerp / (satPara + sunDist));
  680.  
  681.     /* case 2: satellite is on side of the Earth facing away from the Sun */
  682.  
  683.     if (beta >= alpha)    /* satellite can see the Sun */
  684.     {
  685.         if (elev >= ZERO && sunElevation <= TWILIGHT)
  686.             return(VISIBLE);
  687.         else
  688.             return(DAY);
  689.     }
  690.  
  691.     else                  /* satellite is in the Earth's shadow */
  692.         return(NIGHT);
  693. }
  694.  
  695. /******************************************************************************/
  696. /*                                                                            */
  697. /* getDoppler: calculates Doppler shift [Hz]                                  */
  698. /*                                                                            */
  699. /*             the inverted sign on the uplink Doppler shift takes into       */
  700. /*             account the fact that the ground station rather than the       */
  701. /*             satellite needs to adjust the frequency                        */
  702. /*                                                                            */
  703. /******************************************************************************/
  704.  
  705. void getDoppler()
  706.  
  707. {
  708.     if (satElevation > TRACKLIMIT1)
  709.     {
  710.         downlinkDopp = -(downlinkFreq + freqOffset) * rangeRate/CVAC;
  711.         uplinkDopp   =  (uplinkFreq + freqOffset*xponderSign) * rangeRate/CVAC;
  712.     }
  713.  
  714.     else
  715.     {
  716.         downlinkDopp = ZERO;
  717.         uplinkDopp   = ZERO;
  718.     }
  719.     
  720.     return;
  721. }
  722.  
  723. /******************************************************************************/
  724. /*                                                                            */
  725. /* getPathLoss: calculates path loss for given downlink and uplink            */
  726. /*              frequency [Hz], range [km], and velocity of light [km/s]      */
  727. /*                                                                            */
  728. /******************************************************************************/
  729.  
  730. void getPathLoss()
  731.  
  732. {
  733.     if (downlinkFreq > ONEPPM && satElevation > TRACKLIMIT1)
  734.         downlinkLoss = 20.0 * log10(FOURPI * satRange / CVAC * downlinkFreq);
  735.     else
  736.         downlinkLoss = ZERO;
  737.  
  738.     if (uplinkFreq > ONEPPM && satElevation > TRACKLIMIT1)
  739.         uplinkLoss   = 20.0 * log10(FOURPI * satRange / CVAC * uplinkFreq);
  740.     else
  741.         uplinkLoss   = ZERO;
  742.  
  743.     return;
  744. }
  745.  
  746. /******************************************************************************/
  747. /*                                                                            */
  748. /* getSquintAngle: calculates squint angle                                    */
  749. /*                                                                            */
  750. /******************************************************************************/
  751.  
  752. void getSquintAngle()
  753.  
  754. {
  755.     double vectLen;
  756.  
  757.     if (attitudeFlag)
  758.     {
  759.         vectLen     = absol(satSquintGl);
  760.         squintAngle = acos( -(scalar(satAttGl,satSquintGl)) / vectLen);
  761.     }
  762.  
  763.     else
  764.     {
  765.         squintAngle = ZERO;
  766.     }
  767.  
  768.     return;
  769. }
  770.  
  771. /******************************************************************************/
  772. /*                                                                            */
  773. /* getPhase: gets orbital phase of the satellite, i.e. the mean anomaly in    */
  774. /*           units of usually 0 to 256                                        */
  775. /*                                                                            */
  776. /******************************************************************************/
  777.  
  778. void getPhase()
  779.  
  780. {
  781.     satPhase = reduce(meanAnomaly/TWOPI*maxPhase + perigeePhase,ZERO,maxPhase);
  782.     return;
  783. }
  784.  
  785. /******************************************************************************/
  786. /*                                                                            */
  787. /* getMode: gets satellite mode                                               */
  788. /*                                                                            */
  789. /******************************************************************************/
  790.  
  791. void getMode()
  792.  
  793. {
  794.     int curMode;
  795.  
  796.     sprintf(modeString,"  ");
  797.  
  798.     if (numModes > 0)
  799.     {
  800.         for (curMode = 0; curMode < numModes; curMode++)
  801.         {
  802.             if (satPhase >= modes[curMode].minPhase && 
  803.                 satPhase <  modes[curMode].maxPhase)
  804.             {
  805.                 sprintf(modeString,"%s",modes[curMode].modeStr);
  806.             }
  807.         }
  808.     }
  809.  
  810.     return;
  811. }
  812.  
  813. /******************************************************************************/
  814. /*                                                                            */
  815. /* getShuttleOrbit: gets orbit number for the space shuttle                   */
  816. /*                                                                            */
  817. /* convention for element sets from NASA  : add -1                            */
  818. /*                             from NORAD : add  0                            */
  819. /*                                                                            */
  820. /******************************************************************************/
  821.  
  822. void getShuttleOrbit()
  823.  
  824. {
  825.     if (satTypeFlag != STS)
  826.     {
  827.         stsOrbit      = 0.0;
  828.         stsOrbitNum   = 0L;
  829.         stsOrbitFract = 0.0;
  830.         return;
  831.     }
  832.  
  833.     stsOrbit      = curOrbit + curArgPerigee * CRREV;
  834.     stsOrbit     += (!strcmp(elementType,"NASA")) ? -1.0 : 0.0;
  835.     stsOrbitNum   = (long) stsOrbit;
  836.     stsOrbitFract = 100.0 * modf(stsOrbit,&dummyd);
  837.  
  838.     return;
  839. }
  840.  
  841. /******************************************************************************/
  842. /*                                                                            */
  843. /* End of function block satcalc.c                                            */
  844. /*                                                                            */
  845. /******************************************************************************/
  846.